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This paper introduces an adaptive time splitting technique for the solution of stiff evo- 
lutionary PDEs that guarantees an effective error control of the simulation, independent 
of the fastest physical time scale for highly unsteady problems. The strategy considers 
a second order Strang method and another lower order embedded splitting scheme that 
takes into account potential loss of order due to the stiffness featured by time-space 
multi-scale phenomena. The scheme is then built upon a precise numerical analysis of 
the method and a complementary numerical procedure, conceived to overcome classical 
restrictions of adaptive time stepping schemes based on lower order embedded methods, 
whenever asymptotic estimates fail to predict the dynamics of the problem. The perfor- 
mance of the method in terms of control of integration errors is evaluated by numerical 
simulations of stiff propagating waves coming from nonlinear chemical dynamics models 
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as well as highly multi-scale nanosecond repetitively pulsed gas discharges, which allow 
to illustrate the method capabilities to consistently describe a broad spectrum of time 
scales and different physical scenarios for consecutive discharge/post-discharge phases. 

Keywords: Time adaptive integration; error control; operator splitting; reaction-diffusion; 
multi-scale reaction waves; multi-scale discharge. 

AMS Subject Classification: 65G20, 65M15, 65Z05, 65L04, 35K57, 35A35, 35C07 

Dedication 

Cet article est dedie a la memoire de Michelle Schatzman. Specialiste des methodes 
de decomposition d'operateur, sa grande clairvoyance scientifique lui a permis 
d'orienter plusieurs chercheurs debutants sur ce sujet a un moment ou il pouvait 
sembler acheve. Michelle aimait dire qu'il n'y a pas de frontiere entre les branches 
des mathematiques et que seule une grande culture permet de naviguer dans cette 
foret et d'y trouver les bonnes techniques pour resoudre un probleme. Ce travail est 
un hommage; a la croisee des mathematiques et de leur applications effectives, il 
tente d'illustrer cette assertion. Michelle, ton dynamisme, ton humour et ton plaisir 
a parler mathematiques nous manquent. 

1. Introduction 

Numerical simulations of multi-scale phenomena are commonly used for modeling 
purposes in many applications such as combustion, plasma discharges, chemical va- 
por deposition or air pollution modeling. In general, all these models raise several 
difficulties created by the high number of unknowns, the wide range of temporal 
scales due to large and detailed chemical kinetic mechanisms, as well as steep spatial 
gradients associated with localized fronts of high chemical activity. In this context, 
faced with the induced stiffness of these time dependent problems, a high performing 
numerical strategy for multidimensional simulations considers a time operator split- 
ting with dedicated high order time integration methods for reaction and diffusion 
problems, in order to exploit efficiently the special features of each problem. Such 
a numerical strategy for time discretization has been presented in [9] and extended 
in [8] with multircsolution techniques for adaptive space discretization. The main 
idea is to use a second order Strang scheme to solve independently reaction and 
diffusion problems in three successive fractional steps, taking into account that for 
multi-scale phenomena better performances are usually expected while ending the 
splitting scheme by the part involving the fastest scales, as it has been proven in [5]. 
Therefore, based on these theoretical results and on the construction of the splitting 
solver, this strategy provides an accurate resolution of such stiff problems even for 
splitting time steps much larger than either the fastest time scales involved in the 
source terms or the time step restrictions related to spatial grid discretizations. 

Up to our days, fixed splitting time step schemes have been largely used in the 
literature [161 241 22] . and the relevance of our numerical strategy [91 8] has been 
evaluated in the framework of stiff reaction waves for which a constant splitting 
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time step is more than reasonable to precisely describe the global coupling of the 
split phenomena. However, such a fixed time stepping strategy would surely lead 
to major difficulties and limitations for problems describing highly non stationary 
models with very different dynamics such as flame ignition and propagation or 
repetitively pulsed plasmas discharges [23] , all the more in the framework of large 
scale simulations. It is thus essential to be able to dynamically adapt splitting 
time steps for the simulation of such multi-scale problems with strongly evolving 
dynamics. 

In order to guarantee a precise description of the coupled multi-scale phe- 
nomenon, this splitting time step adaptation strategy must rely on a local error 
estimate, which can be obtained by considering a lower order embedded method. 
This is a common practice for ODEs numerical solution [13] , which yields very effi- 
cient and eventually high order methods for which time steps can dynamically adapt 
according to a given tolerance, to sufficiently small values in order to cope with the 
fastest time scales of the problem. However, it is well known that for stiff problems 
and larger accuracy tolerances, the order of the methods can degenerate, yielding 
non reliable error estimates and possibly, much larger global errors than expected 
by the given tolerance. Such a scenario will be all the more valid in the framework of 
the resolution of PDEs where fine grid and large gradients coupled with stiff source 
terms lead to especially stiff problems. In particular, our numerical strategy |91 8j is 
built in such a way that the main source of error is the splitting error, each building 
block relying on high order adaptive and dedicated numerical methods; therefore, 
it is essential not only to construct a reliable splitting error estimate, but also to 
guarantee an effective error control within the so claimed accuracy tolerance. 

In this article, we present a novel strategy to control the local splitting error 
with two different splitting schemes, the first one is a second order Strang technique 
whereas the second one considers a shifted Strang formula, built with a e-shift in 
time of the classical Strang formula. This second method is embedded because the 
first substep is common to both methods to reduce computational cost, and inherits 
from the Strang scheme, stability properties and the same numerical behavior in the 
context of stiff problems; nevertheless, it is only of order one due to the slightly lack 
of symmetry. In the first part of the paper, we conduct a complete error estimate 
of this new splitting method in order to characterize the local error estimate that 
will be computed out of first and second order splitting resolutions. We define then 
a domain of application of the adaptive method in which the local error estimates 
guarantee an effective error control of the solution according to the given tolerance. 
The key issue is related to the evaluation of a maximum splitting time step, called 
the critical splitting time step, as a function of e, for which local error estimates 
are valid. A numerical validation of the theoretical estimates is performed in the 
framework of traveling reaction waves for a simple PDE, for which the threshold and 
critical time steps can be also theoretically estimated and compared with numerical 
results. 

However, in order to extend the numerical strategy to more realistic configu- 
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rations, for which theoretical evaluation of critical time steps is out of reach, we 
develop a complementary and general numerical procedure based on numerical es- 
timates, that allows to establish the domain of application of the method by simul- 
taneously choosing the appropriate e for a given tolerance. This procedure is tested 
in the framework of nonlinear chemical dynamics of Belousov-Zhabotinsky (BZ) 
reactions in a very stiff case in both time and space, yielding satisfactory results. 
As a consequence, a final numerical strategy is conceived that considers adaptive 
splitting time steps and that evaluates simultaneously critical time steps as well as 
best-suited e, in order to guarantee error control for a given accuracy tolerance of 
the simulation with splitting time steps as large as possible. The relevance of the 
proposed strategy is first evaluated for the BZ reaction-diffusion equations, whereas 
a more complex problem issued from the simulation of multi-pulsed gas discharges 
involving several dynamics with very different typical time scales, constitutes the 
second test-case. It is shown that for this second very stiff reaction-diffusion sys- 
tem, splitting time steps can cover a range of three orders of magnitude and always 
guarantee a proper respect of the prescribed tolerance. 

The paper is organized as follows: section [5] describes the adaptive time splitting 
strategy; in section [3J we perform the numerical analysis of the proposed method 
and identify the limit of validity of the local error estimate which is at the heart 
of the adapting procedure. Section 0] is devoted to the validation of the previous 
theoretical estimates and to a theoretical/numerical study of the critical splitting 
time steps in the context of a ID reaction-diffusion problem featuring traveling 
wave solutions. In section [5] we present the final numerical strategy that includes an 
additional numerical procedure to evaluate critical time steps and suitable e. The 
potential of the method is illustrated for the proposed two test-cases in section [6] 
We end in the last part with some concluding remarks. 

2. Adaptive Time Splitting Method 

Let us first set the general mathematical framework of this work. A class of multi- 
scale phenomena can be modeled by general reaction-diffusion systems of type: 



where f : R m -> R m and u : R X E d -> M m , with a tensor of order dxdxmas 
diffusion matrix D(u). 

In the following we will focus on the simplified case of linear diagonal diffusion, 
for which the elements of the diffusion matrix are written as ■Dtii 2 i 3 (u) = D^di^ 
for some positive indices i\, 1%, so that the diffusion operator reduces to the heat 
operator with some scalar diffusion coefficient Di 3 for component Ui 3 of u. A scalar 
one-dimensional model is considered in order to simplify the presentation, taking 



d t u-d x (D(u)d x ii) = f(u), 
u(0,x) = u (x), 




(2.1) 
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into account that extension into higher dimensions of x or u is straightforward: 

d t u-dlu = f(u), i£l, t > 0, 

u(0, x) — uq(x) 1 i£l, t = 0, 

where / and uq are smooth functions. We denote by T 1 uq the solution of (|2.2p . 

Introducing standard decoupling of the diffusion and reaction parts of (|2.2[) . we 
denote by X 1 uq the solution of the diffusion equation: 

d t u D -d 2 u D = 0, x£R,t>0, (2.3) 

with initial data wd(0, •) = uo(') after some time t; and by Y t ua, the solution of 
the reaction part where spatial coordinate x can be considered as a parameter: 

d t u R = f{u R ), xeR,t>0, (2.4) 

with u R (0, ■) = u (-). 

The two Lie approximation formulae of the solution of system (|2.2j) are then 
defined by 

L\u = XVuo, L\u Q = Y t X t u , (2.5) 
whereas the two Strang approximation formulae [25 26 are given by 

S\u Q = X t/2 y*X t/2 u , Slu = Y t/2 X t Y t/2 u . (2.6) 

It is well known that Lie formulae (|2.5[) (resp. Strang formulae (12.61) ) are an ap- 
proximation of order 1 (resp. 2) of the exact solution of (I2.2p . Higher order splitting 
schemes are also possible. Nevertheless, the order conditions for such composition 
methods state that either negative time substeps or complex coefficients or non con- 
vex combinations are necessary [13] . The formers imply usually important stability 
restrictions and more sophisticated numerical implementations. In the particular 
case of negative time steps, they are completely undesirable for PDEs that are 
ill-posed for negative time progression. 

An adaptive time stepping strategy is based on a local error estimate which can 
be obtained by using two schemes of different order, in this case S[ or S^, locally 
of order 3, and L\ or L\ : locally of order 2. For instance, the Embedded Split-Step 
Formulae given in [17] consider S\ and L\ or S\ and L\, noticing that 

where Y 1 I 2 uq is also used to compute S\ uq. Nevertheless, in the context of multi- 
scale phenomena, order reductions may appear due to short-life transients associ- 
ated with the fastest variables when one considers splitting time steps larger than 
the fastest scales. It has been proved in [5] that better performances are expected 
while ending the splitting scheme by the part involving the fastest time scales of 
the phenomenon. In particular, in the case of linear diagonal diffusion problems, 
no order loss is expected for the L\ and S\ schemes when fast scales are present in 
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the reactive term. Therefore, the embedding procedure must be carefully conceived 
taking into consideration these theoretical studies. 
We introduce a shifted Strang formula 

Sl e u = Y^-^X'Y^+^uo, (2.7) 

locally of order 2, due to the lack of symmetry, for e in [—1/2, 0) U (0, 1/2]. In this 
way, a local error estimate is computed based on two solutions for which orders are 
guaranteed and a potential loss of order is simultaneous, following 

fS At u \ ( Y At ' 2 X At Y At / 2 u \ 



S At e u J I Y (1 / 2 -^ At X At Y (V2+ £ ) At U() 



(2- 



for some splitting time step At > 0. Embedding is accomplished as long as e is 
different from —1/2, that is S A *uo different from L\uq. On the other hand, if e is 
equal to 1/2, S A *uo is defined as L\uo, which it is not suitable for stiff configurations 
as it was previously discussed [5]. Therefore, e should be contained in (—1/2,0) U 
(0, 1/2). Shifted S A luo is defined in a similar way and depending on the multi-scale 
character of the problem, it might be the appropriate choice along with S At UQ. 
Taking into account that 

5 2 A *u - S At e u = S At u - T At u + T At u - S A *u , 

= 0(At 3 )+0(At 2 )^0(At 2 ), (2.9) 

for a given accuracy tolerance 77, 

\\S At u - S^uoW <n (2.10) 

must be verified in order to accept current computation with At, while new time 
step is calculated by 



At— = vAt V (2.11) 

y p 2 "o - S^uo\\ 

with security factor < v < 1 close to one. This comes from a classical adaptive time 
stepping procedure for stiff ODEs solution, for which more sophisticated formulae 
than (|2.11|) can be also considered, see [14] for example. 

The error control of these adaptive methods is fully guaranteed as long as the 
orders of both, the main and the embedded integration methods, remains valid. This 
is the case for small enough time steps for which asymptotic theoretical estimates 
hold, but remains an open problem for larger time steps for which the validity of the 
formers is assumed. This is a key point in this work, because we propose not only a 
new splitting strategy with adaptive time steps as described in this section, but we 
aim also at applications for which splitting time steps may go beyond the fastest 
scales associated with each subproblem in order to obtain important computational 
savings. Therefore, a technique that guarantees consistently error control for all 
possible separation scales must be pursued, but first of all, a detailed numerical 
analysis of the method must be performed. This is the goal of the following part. 



1 
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3. Numerical Analysis of the Method 

In this part, we develop the numerical analysis of the proposed method. It is mainly 
based on the theoretical study of the introduced shifted Strang formula (|2.7[) and 
the domain of validity of the local error estimates. However, first of all, we introduce 
the Lie formalism which will be used as mathematical tool of analysis. 



3.1. The Lie operator formalism 

We introduce the Lie operator formalism in order to generalize the exponential of 
a linear operator in the context of nonlinear operators. Let A be a Banach space, 
To > and F, an unbounded nonlinear operator from D(F) C X to X, we consider 
the general autonomous equation: 




(3.1) 



u'(t) = F(u(t)), 
u(0) = u , 

The exact solution of this evolutionary equation is (formally) given by 

u(t) = T 1 uq, 0<t<T o (3.2) 

where T* is the semiflow associated with (|3.1|> : in particular we can set F(u) — 
d^.u + f(u) as in (|2.2p . The Lie operator Dp associated with F is then a linear 
operator acting on the space of operators defined in X [131 6j . More precisely, for 
any unbounded nonlinear operator G from D(G) C X to X with Frechet derivative 
G', Dp maps G into a new operator DpG, such that for any v m X: 

(D F G)(v) = G'(v)F(v). (3.3) 

Hence, by induction on n with solution u of (|3.1[) . we obtain 

^G(u(t)) = (D F G)(u(t)), 
and a formal Taylor expansion yields 

G(u(t)) = Y^(^G(u(t)) 



71 = 



t=o 




(e tD -G)u . (3.4) 



If we now assume that G is the identity operator Id, we obtain 

u(t) = T l uo = (e tDF ld) u a . 

Therefore, the Lie operator is indeed a way to write the solution of a nonlinear 
equation in terms of a linear but differential operator. Following (|3.4p . an important 
result obtained by Grobner in 1960 [12], considers the composition of two semifiows 
T\ and T| associated with F\ and F^ for any d in X: 

T(Tiv = (e sD ^Tl) v = (e sDF ^e tD ^ld) v. 



i 



January 19, 2013 23:50 WSPC/INSTRUCTION FILE AdaptiveSplit- 
ting' revised 



8 Descombes, Duarte, Dumont, Louvet, Massot 

3.2. Error analysis 

In this paragraph, we conduct the error analysis of the approximation of T* by 
5| e in a linear framework. Then, we extend these results to a general nonlinear 
configuration given by problem (|2.2j) . using the Lie operator formalism. General 
estimates for the approximation of T* by are a l so drawn. We end in the last part 
with a mathematical study that shows the domain of application of the method 
described in §[2l for which an effective error control is guaranteed within an accuracy 
tolerance. To simplify the notations in what follows, we will denote S\ by S l and 
S{ e by 

Assume that A and B are linear bounded operators and define 

Slu = e (l/2-e)« e tB e (l/2 +e) a Uo 

as an approximation of e '( A + s ). The following theorem gives the expansion in pow- 
ers of t of the difference between e l ( A + B ) and St. We recall the definition of the 
brackets between A and B: [A, B] = AB - BA. 

Theorem 3.1. Assume that A and B are linear bounded operators, for t and e 
small enough, the following asymptotic holds 

e^+^uo-Stuo = -et 2 [A, B]u + ^([A, [A, B}] +2[B, [A, B}] )u + O(et 3 ) + O(t 4 ). 

Proof. Proof is straightforward by using the Taylor formula with integral remain- 
der for a linear bounded operator A: 



tA . t 2 A 2 t 3 A 3 [* (t-s) 3 td . 

e tA = Id + tA+ — + — + / v . ; A 4 e sA ds. 



□ 



6 Jo 6 

We extend now the previous theorem to our nonlinear framework given by (|2.2[) . 
In order to do this, we introduce the spaces C°°(K) of functions of class C°° on R, 
and C£°(K) of functions of class C°° on R and bounded over R. We consider also 
the Schwartz space <S(R) defined by 

5(R) = {g e C*°°(R) | sixp\v ai d* 2 g(v)\ < oo for all integers a u a 2 }; 

and we define the space Si (R) , made out of functions v belonging to (R) such 
that v' belongs to 5(R). Let us consider now equation (|2.2p and give the expansion 
in powers of t of the difference between T* and 5 1 *, given by (|2.7[) . 



Theorem 3.2. Assume that uq belongs to S\ (R) and that f belongs to C c 
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For t and e small enough, the following asymptotic holds 

T t u -S t e u = -et 2 f"(u ) (^f 2 

+^(/'M/"M + /M/ (3) M) (jlf)' 

12" 

-¥ //,(uo) (S?) 2+0(rf3)+0( * 4) ' (3 - 5) 

Proof. We introduce the two Lie operators Da and Df associated with d 2 and / 
and write 

T'uo - Slu = (e'^+^Id) u - ( e (i/24*)*D /e *D Ae (i/2- e )tD /Id ^ UQ _ 

With Theorem l3.ll we can deduce that 

T'uo - Stuo = -et 2 ([Df, Da] Id) u + ^ ([D f , [D f ,D A ]]ld) u 
t 3 

+^ ([DAADf^A^uo + Oiet^ + Oit 4 ). (3.6) 

We are not interested in giving the exact form of the terms 0(et 3 ) and 0(t A ), but 
these terms can be computed following the same technique developed in [6] . For the 
term in 0(t 2 ), we have by definition and with (|3.3|) . 

{[Df> Da] Id) u Q = (D f (D A ld) - D A (D / Id)) u , 

d 2 u 



(D A id)'K)/K) - (D f idy(u y 



dx 2 



The last term is by definition the Lie bracket between d 2 and /, a simple compu- 
tation shows that 



d 2 /M s d 2 u „ ,(du \ 2 , (Puq , d 2 u 



Furthermore, 



([D f , [Df, Da]] Id) (uo) = (f(u )f"(u ) + /( Uo )/ (3) M) (^j 
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and 



dx J V dx ) dx 2 

8 2 uq 



([D A , [D f ,D A] ]U) u = -/ (4) K) ('^V - 4 / (3) M ( d ~^ V d2U ° 

-2/"(-o) , dx2 

All the terms are now computed and this concludes the proof of Theorem 13.21 □ 

For e = 0, the next corollary follows directly. 

Corollary 3.1. Assume that Uq belongs to S\ (R) and that f belongs to C°°(R). 
For t small enough, the following asymptotic holds 

T'uo - fl* Uo = ^(/'(uo)/"(«o) + JM/ (3) M) (^j 



8uq\ 2 9 2 wo 
dx 2 



6 \ dx 2 

From (|3.5[) and (|3.7|) . we can see that 



£/>o)(^] +0(t% ,3.7) 



5*u - S*t*o - rf 2 /"K) ( ^ J +0(et 3 ), (3.8) 



and thus, 



T*w ~ 5*m = T'up - g'ug + S"u - S^uo . (3.9) 

0(* 3 ) 0(et 2 ) 

Therefore, we are sure that the real local error of the method, T*uo — S^tto, will be 
bounded by the local error estimate, err = S t u — S^Wq, when for a given e, 

T*m - Sluo « C(t 2 ) (3.10) 

is verified into Q3.9P ; that is, when the embedded method is really of lower order as it 
was assumed in (|2.9[) . This will be always verified for small enough time steps t, for 
which T*u — 5*Mo ~ 0(t 3 ) < err ss 0(et 2 ) is guaranteed. Nevertheless, for larger 
time steps, err will fail to properly predict T'mo — S Uo since we will eventually 
have T'm - SPuo ~ 0(t 3 ) > err w £>(ei 2 ). When this happens, (|3.10|l is no longer 
true and the previous estimates show that we will rather have T*«o — S^uq ~ 0(t 3 ), 
and assumption (|2.9p will no longer hold. 

In order to overcome this difficulty, we must therefore estimate a critical time 
step t* > such that for all t in [0, t*], (|3.10l) is guaranteed for a given e. This will 
imply that Strang local error, T'uo — S*uq, will be indeed bounded by the local error 
estimate, err, and that an effective error control will be achieved for err smaller 
than a given accuracy tolerance r\. Finally, a suitable choice of e can be also made 
since t* is related to e following (|3.9I) . 
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A natural strategy to predict this critical t* will rely on the previous theoretical 
estimates and on a more precise knowledge of the structure of the solutions of the 
PDEs; this is for instance illustrated in the next part in the context of traveling 
wave solutions. 



4. Application to Reaction Traveling Waves 

In this part, we will confront the previous theoretical study to a simple reaction 
diffusion problem that admits self-similar traveling wave solutions such as the KPP 
equation [IB]. The main advantages of considering this kind of problems are that 
analytic solutions exist and that the featured stiffness can be tuned using a space- 
time scaling. Therefore, it provides a first numerical validation of the numerical 
estimates of the method and an evaluation of its domain of application; and on 
the other hand, a detailed study can be conducted on the impact of the stiffness 
featured by propagating fronts with steep spatial gradients. 

In what follows, we recast previous estimates in the context of these reaction 
traveling waves, to then deduce an estimate of the time step t* that defines the 
limit of application of the method for which local error estimates yield effective 
error control. We end with a numerical validation of the theoretical results in the 
context of the resolution of KPP model. 



4.1. Numerical estimates 

We are interested in the propagation of self-similar waves modeled by parabolic 
PDEs of type: 



d t u -Dd 2 x u = kf(u), t > 0, 

u(0,x) = uq(x), x G K, t = 0, 



(4.1) 



with solution u(x, t) = uq(x — ct), where c is the steady speed of the wavefront, and 
D and k stand respectively for diffusion and reaction coefficients. 

Considering Theorem 13.21 we obtain the following estimate for system (|4.1|) . 



Corollary 4.1. Assume that uq belongs to S\ (M) and that f belongs to C°°(l 
For t and e small enough, the following asymptotic holds 

T*uo - S*«o = ~ekDt 2 f"(u ) " 

+^(/'M/>o) + /M/ (3) M) (^) 2 

kD 2 t 3 .„. . fd 2 u 



(, /M U^' + ° {Et (4 - 2) 



I 
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Proof. Proof follows directly from demonstration of Theorem 13.21 using (|3.6[) and 
considering that 

[D k f,D DA ] = kD[D f ,D A ], 
[[D kf ,D DA ],D DA ] = kD 2 [[D f ,D A ],D A ], 
[[D kf ,D DA ],D kf ] = k 2 D[[D f ,D A ],D f ], 



where Dd A and Df-f are the Lie operators associated with DO 2 and kf. 



□ 



On the other hand, if we now consider system (|4.1|) with k = 1 and D = 1, the 
following corollary establishes t* > such that for all t in [0, t*] (|3.10p is guaranteed 
for a given e. 

Corollary 4.2. Assume that uq belongs to Si (R) and that f belongs to C°°(K). 
For a given e small enough , define 

(duo" 2 



Mi 



/"M 



\ dx 



(4.3) 



L' 2 



and 



M 2 = 



f'(uo)f"(u ) + f(u )fW(u ) fdu^ 2 
24 V 9x 

f^(u ) {dugV&Uo _ />) /9^o 

dx J dx 2 6 \ dx 2 



/ (4) ("o) 
12 V 



define t* by 



t*M 2 = zM x . 

For all t such that < t < t* then 



(4.4) 



(4.5) 



C(i 2 ). 



In a general case, if evaluation of the derivatives of uq and / is feasible, it 
is then possible to predict the domain of application of the method, [0,t*], for 
a given e based on the previous result. In the particular case of traveling wave 
solutions for (14. ip . diffusion and reaction coefficients, D and k, might be seen as 
scaling coefficients in time and space. A dimensionless analysis of a traveling wave, 
as shown in [11], can be then conducted considering a dimensionless time r and a 
dimensionless space r with 

r = kt and r = (k/D) 1/2 x. 

This analysis allows to find a steady velocity of the wavefront, 

c = x t oc (Dk) 1/2 , 

whereas the sharpness of the wave profile is measured by 



^La X «(fc/^) 1/2 - 



(4.6) 
(4.7) 
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Therefore, condition Dk = 1 implies constant velocity for all k = 1/D but greater 
k (or smaller D) implies higher spatial gradients, and thus, stiffcr configurations. 

This study gives complementary information on the solution of (|4.ip and in 
particular, when condition Dk = 1 is satisfied, it allows to deduce from Corollary 
lL2l 

M*M 2 = eM 1 , (4.8) 

with Mi and A'h given by (|4.3|) and (|4.4[) . Therefore, stiffer configurations given by 
the presence of steeper spatial gradients will restrain the application domain of the 
method, according to (|4.8I) . Nevertheless, for larger gradients, smaller time steps are 
also required for a given level of accuracy and hence, we can expect a simultaneous 
reduction of both critical and accurate splitting time steps. 

4.2. Numerical illustration: KPP equation 

Let us recall the Kolmogorov-Petrovskii-Piskunov model. In their original paper 
[18] . these authors introduced a model describing the propagation of a virus and 
the first rigorous analysis of a stable traveling wave solution of a nonlinear reaction- 
diffusion equation |11) . The equation is the following: 

d t u- Dd 2 x u = ku 2 {l-u), (4.9) 

with homogeneous Neumann boundary conditions. We consider a ID discretization 
with 5001 points on a [—70, 70] region for which we have negligible spatial discretiza- 
tion errors with respect to the ones coming from the numerical time integration. 

The description of the dimensionless model and the structure of the exact solu- 
tion can be found in [IT] where the dimensionless analysis shows that in the case 
of D = 1 and k = 1, the velocity of the self-similar traveling wave is c = 1/ ' \[2 and 
the maximal gradient value reaches l/y/32. The key point of this illustration is that 
the velocity of the traveling wave is proportional to (k D) 1 / 2 , whereas the maximal 
gradient is proportional to (k/D) 1 / 2 . Hence, we consider the case kD = 1 for which 
one may obtain steeper gradients for the same speed of propagation. 

Throughout all this paper, exact solution T*uo will be approximated by the res- 
olution of the coupled reaction-diffusion problem performed by the Radau5 method 
[15] with fine tolerances, r\R a dau$ = 10~ lf} . This solution will be referred as the refer- 
ence or quasi-exact solution. Strang approximations S t UQ and SIuq will be computed 
with a splitting technique recently introduced 18 1 9] , which considers Radau5 [15] to 
solve locally point by point the reaction term; and the ROCK4 method [1] for the 
diffusion problem. Radau5 [TS] is a fifth order implicit Runge-Kutta method exhibit- 
ing A- and ^stability properties to efficiently solve stiff systems of ODEs, whereas 
ROCK4 [JJ is formally a fourth order stabilized explicit Runge-Kutta method with 
extended stability domain along the negative real axis, well suited to numerically 
treat mildly stiff elliptic operators. Both methods implement adaptive time step- 
ping techniques to guarantee computations within a prescribed accuracy tolerance. 
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In order to properly discriminate the previously estimated splitting errors from those 
coming from temporal integration of the substeps, we consider also fine tolerances, 

?7R a dau5 = f?ROCK4 = 10~ 10 . 




Fig. 1. KPP equation with k = 1. Local L 2 errors for several splitting time steps At and e = 0.05 
(top left), 0.005 (top right) and 0.0005 (bottom left). Bottom right: critical splitting time steps 
At* obtained when ||T A *tio — S^uqW^ ~ ||S' At Mo — S^uqW^ in the numerical tests. 



Figures Q] and [2] show L 2 errors between T l uo, S'lto and S^uo solutions for k = 1, 
k = 10 and k = 100 respectively, and several s. Notice that estimates (I3.5[) . (|3.7|) 
and ()3.8|) for all three errors in (|3.9p are verified and in particular, for At larger than 
critical At*, the estimated error err = \\S At UQ — S^uollz, 2 is no longer predicting 
the real local error given by T'uo — S'lto. 

With these results, we can also compare real At*, obtained when ||T At u — 
"S^'woIIl 2 ~ \\S Ai uq — S^'uoIIl 2 in the numerical tests, with theoretically esti- 
mated At* following ()4.8|) . Table Q] summarizes these results where computation of 
estimated At* in (|4.8|) is given by the computation of Mi and with Maple® 
according to (|4.3p and (14. 4|) . A really good agreement can be observed even though 
theoretical results underestimate the real values. The loss of order depicted by the 
numerical results, is due to the influence of spatial gradients in the solution, as 
it was proven in [3]. This explains the error of the predicted critical At* in (|4.8[) 
whenever one gets close to the order loss region. 

Numerical results show also that ||S' A 'uo — S^uoIIl 2 oc s according to (|3.8|) 
and consequently, At* oc e; therefore, the working region or domain of application 
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Fig. 2. KPP equation with k = 10 (top) and k = 100 (bottom). Local L 2 errors for several 
splitting time steps At and e = 0.05 (left). Right: critical splitting time steps At* obtained when 
\\T At uo — S At uo||^2 ~ ||>S At Jio — 5 a *uoIIl2 m the numerical tests. 



of the method, At < At*, depends directly on the choice of e as it can be seen 
in Table [TJ Finally, in the context of traveling waves, these numerical experiments 
show that At* cx k^ 1 cx l/||<9iio/<9x|| 00 according to Tabled] hence, application 
domains are reduced for stiffer configurations but numerical results show also that 
smaller time steps are required for the same level of accuracy. These conclusions 
are easily extrapolated to more general self-similar propagating waves. 



Table 1. KPP equation. Comparison between real At* cal , obtained when ||T At uo — 5 At «o|| i 2 ~ 
HS^'uq — S^'uoIIl 2 in the numerical tests, and theoretically estimated At* st following 114.81 1. 





e = 0.05 


£ = 0.005 


e = 0.0005 


k = 


1 


At* 


2.783 


0.1274 


1.17 x 10~ 2 






At* 


1.107 


0.1107 


1.11 x 10~ 2 


k = 


10 


At* 

^real 


0.2803 


1.29 x 10~ 2 


1.19 x 10~ 3 






At* 


0.1107 


1.11 x 10~ 2 


1.11 x 10~ 3 


k = 


100 


At* 

^real 


4.33 x 10~ 2 


2.12 x 10~ 3 


1.92 x 10~ 4 






At* 


1.11 x 10~ 2 


1.11 x 10~ 3 


1.11 x 10~ 4 
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5. Construction of the Numerical Strategy 

We have presented in § [21 a time adaptive numerical scheme fully based on theo- 
retical error estimates developed in § [31 We have also studied the necessary general 
conditions in order to guarantee an effective error control based on local error esti- 
mates. In particular, this has been shown in the case of reaction traveling waves in 
§[H for which theoretical studies give us some insight into the PDE solution. Never- 
theless, this is not always possible and it is usually difficult to carry out such kind 
of analysis for more realistic models. Therefore, based on the theoretical analysis 
and previous illustrations on the influence of the various parameters of the scheme, 
a general numerical procedure that completes the adaptive scheme defined in § [21 
is introduced in the following. 

In a first part, we will settle the theoretical framework and the numerical proce- 
dure needed to estimate t*, and to define the appropriate e. This will be illustrated 
by numerical tests performed on a more complex model of time-space stiff propa- 
gating waves. These theoretical and numerical studies will allow to define, at the 
end, a final numerical strategy. 

5.1. Numerical procedure to estimate critical t* and e 

Let us consider general system (|2.2|) , based on theoretical estimates (|3 . T[) and (|3.8[) , 
we can write 

S At u - T At u = C At 3 , (5.1) 

where C = d{u ) + 0(At 4 ), and 

S At u - S At u Q = eC £ At 2 , (5.2) 

where C E = Ca(uo) + 0(e, At 3 ); the dependence of C e on e is only given in the 
higher order terms and it is thus neglected. 

For a given e, in the same spirit as Corollary 14. 2\ we search for a critical At* 
such that 

\\S At u - T At u \\ < \\S At u Q - S At u \\ (5.3) 
for all At < At* . According to (|5.1|) and (|5.2I) . we have then the following estimate: 

At* » f£. (5.4) 

For a given e, this gives an upper bound for the time steps for which the local error 
estimate, err — \\S Uq — S At uo\\, is properly estimating the real Strang local error, 
\\S At u - T A *w ||, following 

In particular, when At — > At*, we have that err \\S At u ~T At u \\, and the 
local error estimate is predicting more accurately the real error of integration. The 
critical time step, At*, is directly related to e through (|5.4p as we have already 
shown in the previous numerical results in § 14.21 Therefore, a suitable e will define 
a critical At* such that the estimated splitting time steps At for a given tolerance r\ 
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will be close enough to critical At* , in order to avoid an excessive overestimation of 
the Strang local error and thus, larger time steps can be chosen for a given accuracy 
tolerance rj. 

In order to compute At* for a given e, we must first estimate Co in (|5.4j) . 
since C £ is computed out of the local error estimate, err, for known At and e in 
(|5.2j) . Estimating Co amounts to directly estimate Strang local error through (|5.1[) 
and thus, the accuracy of the simulation might be controlled in this way without 
relying on a local error estimate as proposed in the embedded method strategy in 
§ [2j Nevertheless, as we will see in the following, in order to estimate Co and the 
Strang local error, we must define new local estimators and a numerical procedure 
that becomes rapidly very expensive if we want to implement such error control 
technique. Therefore, we must rely on a local error estimate given by a less expensive 
strategy for which the computation of Co is only performed from time to time to 
guarantee the validity of local error estimates. 

The next Lemma will be useful to define the numerical procedure to estimate 

C . 

Lemma 5.1. Let us consider system 12. 2\) and assume a local Lipschitz condition 
for f: 

ll/(«)-/(w)ll<A||u-i;||. (5-5) 
For a finite At the following holds 

\\T At u ~T At v \\ <w||«o-«b||, (5-6) 
with u) = 1 + nAt for small enough At. 

Proof. Using Duhamel's formula for (|2.2[) yields 

T'uo - T*«o = e td * (u ~ v ) + f e ( '" s ^ (f(T s u ) - f(T s v )) ds. (5.7) 

Jo 

Taking norms and applying recursively (|5.7|) . 

||r'uo - T^woll < || wo - "o|| + A / ||T s it - T s v \\ ds, 

Jo 

< e xt \\u -v Q \\, (5.8) 
proves (|5.6p for t — At finite. □ 

If we define a local estimator, ei = S aiAt u - S blAt (S clAt u ), such that ai = 
b\ + Ci, we obtain that 

S blAt (S clAt u Q ) - T aiAt u = S blAt (S clAt u ) - T blAt {S clAt u ) 

+T blAt (S ClAt u Q ) - T blAt (T ClAt u ), 

= Csci At u °lAt 3 

+T blAt {S ClAt u ) - T blAt {T ClAt u ), (5.9) 
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where C s ^t Uo = Ci(5 ClAt u ) + 0{At 4 ). Therefore, assuming that Cgc 1 &t UQ m Cq 
and considering Lemma [5.1[ it follows from the difference between (|5.ip at a\At 
and (|53t : 

\\ ei - (a 3 , - bl)C At 3 \\ < uj\\T^ At uo - S^ At u \\, 

< ujC clAt 3 . (5.10) 

Hence, defining a second local estimator, e 2 = S a2At uo — S b2At (S C2At uo), such 
that 02 = 62 + C2, we obtain a second expression similar to (|5.10l) with ei and 
(02,62,02), and we can estimate Co and u. In particular, we notice that b\ should 
be close to 62 in order to better approximate to into (|5.6p and (|5.10p . and that c\ 
and C2 should also be small enough to guarantee Cs ol At U0 w Co and Cgc 2 At UQ m Cq. 
On the other hand, to optimize the required number of extra computations from 
a practical point of view, we can use the estimator e2 to compute estimator e\ by 
setting ei2 = Ci, and we can also fix a\ = 1 so we can use S aiAt uo for the time 
integration of the problem. In this way, the extra computations needed to compute 
local estimators ei and e 2 will be given by S c * At u , S b2At (S C2At u ), S ClAt u and 
S blAt (S ClAt uo) within a time step At. Then, we will be able to compute oj and Co, 
by solving two expressions of type (|5.10[) . The next numerical example illustrates 
the validity of this numerical procedure. 

5.2. Numerical example of evaluation of critical t* : BZ equation 

We are concerned with the numerical approximation of a model of the Belousov- 
Zhabotinski reaction, a catalyzed oxidation of an organic species by acid bromated 
ion (for more details and illustrations, see [TU])- We thus consider the model intro- 
duced in [TT] and coming from the classic work of Field, Koros and Noyes (FKN) 
(1972), which takes into account three species: HBrC>2 (hypobromous acid), bromide 
ions Br~ and cerium(IV). Denoting by a = [Ce(IV)], b = [HBrC^] and c = [Br - ], 
we obtain a very stiff system of three partial differential equations: 

d t a - D a &la 

d t b-D b d 2 x b 
d t c - D c dlc 

with diffusion coefficients D a , D b and D c , and some real positive parameters /, 
small q, and small e, a, such that fj, <C e. 

The dynamical system associated with this system models reactive excitable 
media with a large time scale spectrum (see [TT] for more details). Moreover, the 
spatial configuration with addition of diffusion generates propagating wavefronts 
with steep spatial gradients. Hence, this model presents all the difficulties associ- 
ated with a stiff time-space multi-scale configuration. The advantages of applying 
a splitting strategy to these models have already been studied and presented in jjj . 



-(-qa - ab + fc), 



-(go -06 + 6(1-6)), * 



(5.11) 



= b — c, 
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We consider the ID application of problem (|5.1ip with homogeneous Neumann 
boundary conditions in a space region of [0, 80] with a spatial discretization of 4001 
points, good enough to prevent important spatial discretization errors, and the 
following parameters, taken from [Tl : e = 10~ 2 , ft = 10~ 5 , / = 3 and q = 2 x 10~ 4 , 
with diffusion coefficients D a = 1, Db = 1 and D c — 0.6. Reference solution and 
Strang approximations are defined in the same way as in the KPP application with 
the same tolerances for the time integration solvers. 

First of all, we validate theoretical order estimates (|3.5|) . (|3.7|) and (|3.8p and ver- 
ify relation (|3.9[) . Figure [3] shows L 2 errors between T*-u , <S'«o and Slug solutions 
for several e and the real At;* such that ||T At u - S At u \\ L ^ rj \\S At u - S At u \\ L ^ 
obtained after treating the numerical results. Maximum L 2 error considers the max- 
imum value between normalized local errors for a, b and c variables; in these nu- 
merical tests, it corresponds usually to variable b. 





- s- v „ 

- <Pu 












ii -V "n 
-S?u 


■ 


























V 










- 










-f^- 








f=0 


05 



























s\ 


ii "ii 
- S A 'u 












r v , 


ii -V "ii 

- S^u 
ii V "ii 


• 






















































r=0 


005 



























V 7 -i n" 6 1 rr 5 



10" D 10"° 10"" 10" 3 10" 2 10" 1 



- A 


- s- v „ 
(1 "(I 

- S A 'u 












ii -V H « 

- 

(1 -V "(I 


■ 














































- 








r=0 


)005 





























-=0.0# 










t 
r= 


-0 005 
3 0005 






































1 A, '-° 


0001494 










00001583 




















/ Ai =0.00000151 


i 









10"° 10"° 1 W* 10"" 10" 



Fig. 3. BZ equation. Maximum local L 2 errors for several splitting time steps At and e = 0.05 
(top left), 0.005 (top right) and 0.0005 (bottom left). Bottom right: critical splitting time steps 



At* obtained when ||T At u - S At u \\ L 2 « \\S l 



"Uo 



S At uo\\ L 2 in the numerical tests. 



Let us now define the two sets (ai,&i,Ci) and (02,62,02), and compute local 
estimators e\ and in order to obtain Co according to (|5.10p with At; = At;o = 
10 -5 ; that is a time step for which there is no order loss yet, as seen in Figure 
[31 As it was previously detailed, we consider a\ = 1 and 02 = Ci to avoid some 
extra computations. Furthermore, 62 should be close to b\, and c\ and C2 small 
enough. Setting b\ larger than 1/2 would yield more different 62 since c\ — a^. On 
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the other hand, for b\ smaller than 1/2 we can even set bi = b\ but in this case c\ 
will be larger than 1/2. Therefore, we reach a compromise by setting b\ = 1/2 that 
yields ci — a 2 = 1/2, so we can choose for instance 62 = 2/5 close to b\, and thus, 
c 2 = 1/10. 

With the local error estimate, err = \\S At uo — S^uqWl 2 , for the various time 
steps and several e shown in Figure 02 Figure 2] presents the estimated critical At* 
calculated with (|5.4I) from the estimated Co(Ato) and err. These critical time steps. 
At*, estimated with (|5.4|) are in good agreement with numerically measured At* in 
Figure [31 and depend on the value of e. Hence, the domain of application or working 
region of the method, At < At*, might be settled depending on the desired level 
of accuracy by means of an appropriate choice of e. For instance, if we consider 
the case e = 0.05 in Figure [3j for At = 10~ 6 , the local error estimate is given by 
err « 10 -10 whereas the real Strang local error is ~ 10~ 12 . This overestimation of 
the local error will certainly imply an underestimation in the required size of the 
time steps for a given tolerance. Therefore, for a given tolerance rj a more suitable 
configuration should consider an e such that At « At* in order to reduce excessive 
overestimations of local errors. 
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Fig. 4. BZ equation. Working region of the method At < At* with At* calculated with Co 
estimated at Aio = 10 — 5 and err obtained for several splitting time steps At and e (left). Right: 
predicted Strang error calculated with Co estimated at Ato = 10 _ 5 and locally at several splitting 
time steps At. 



In the illustration shown in Figure |4j Co was estimated in the third order region 
of the method and therefore, all values are well approximated as long as At remains 
in this region. In particular, critical At* will be progressively underestimated for 
larger e and consequently, it will impose smaller time steps for a given tolerance; 
this is already the case for e = 0.05, for which At* is in the transition zone towards 
the lower order region. Even though the computation of Co with small time steps 
will be less expensive, a much more accurate procedure considers current time step 
as shown in Figured) In particular, by estimating locally Co, we are estimating real 
Strang error and thus, At < At* guarantees prescribed accuracy even if asymptotic 
order estimates are no longer verified. This allows to properly extend the domain 
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of application over the whole range of possible time steps for a given accuracy; an 
extremely important issue for real applications for which splitting time steps may 
go far beyond asymptotic behavior including the potential order reduction region 
associated with the stiffness of the problem. 

5.3. Numerical strategy 

Previous studies conducted in § l5.1l and § l5.2l allow to properly complete the adaptive 
splitting strategy introduced in § [2j In this part we conduct the final description of 
the numerical strategy. 

Let us consider general problem (|2.ip for u £ W m , for which we use S\ in (12.61) as 
resolution scheme. Depending on the problem, the adaptive method will be applied 
considering time evolution of I < m variables: u £ M 1 . Let us denote fi; the set of 
indices of these variables. In order to consider only I < m variables, the formers 
must be decoupled of the remaining m — I variables in the reactive term f(u) in 
(|2.1j) . To simplify the presentation, we will only consider e £ (0,£ max ), e max < 1/2. 

We set the accuracy tolerance rj, an initial time step At and initial So, and 
perform the time integration of (12.11) with the Strang scheme S\ and the embedded 
shifted one S*| £ given by (|2 . 7[) . We compute local error estimate err and new time 
step At new according to (|2 . 1 1 [) . If err is smaller than 77, current time step solution 
is accepted and simulation time evolves; otherwise, current solution is rejected and 
the time integration is recomputed with At ncw . In particular, it is better to choose 
rather small At to avoid initial rejections. 

In order to guarantee an effective error control, we define the working region 
At < At* by estimating the corresponding At* for current e. This is done for the 
first time step At and then periodically after N accepted time steps depending on 
the problem, based on the numerical procedure introduced in § 15.11 Computation 
of critical At* is also performed with u, and a rather large initial Eq is suitable to 
initially guarantee At < At*. 

We define then a suitable working region At £ [f3 At* , -f At*] with < j3 < 7 < 1, 
for which splitting time steps are close to At*. A new e is then computed if At is 
much lower than At* (At < /3At*) in order to avoid unnecessary small time steps; 
or if At is very close or possibly larger than At* (At > 7 At*) with 7 close to one, in 
order to increase upper bound of the domain of application. This guarantees that e 
is dynamically computed and properly adapted to the dynamics of the phenomenon. 

Finally, the numerical resolution strategy can be summarized as follows, where 
U £ R mxn stands for the spatial discretization of u over n points, U := (u^'V) 
such that j £ [1, m] and k £ [l,n]. 

• Input parameters. Define accuracy tolerance 77, time domain of study 
[to,T], initial time step At , initial and period of computation of At*: 
N. 

• Initialization. Set iteration counter i = and t = to, U = Uo, At = At , 
e = £q. We define a flag estimate initialized as .false.. Throughout the 
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whole computation, we need to store U, an array of size m x n. 
• Time evolution. If t < T: 



(1) Only if ± 



or estimate is .true . : 



I 

N_ 

Computation of critical At* I: For the sets (ai,&i,ci) and 
(02, 62, C2) with ax = 1 and 02 = Cx, we compute successively: 

- Ui = 5 C2At U , where U is built out of U, U = (u^) jeQl ; 

- Ui = 5 b2At U i; 

- U 2 = 5 ClAt U ; 

- e x = maxjen, ||m 2 : ''' ) - &i ||; 

- U a = 5 blAt U 2 ; 

- estimate is set to . true . . 

These operations needs to store Ux and U 2 , two arrays of size I x n. 

(2) Time integration over At: We compute successively: 

- for each fc G [l,n], it&S = y At /V^); 

- for each fc G [l,n], = F £A * ui£} ; 

- U* = A A 'U*, with U* = f (U new ,Ui); 

- for each fc e [l,n], - rCVa-^At^fc). 

- for each fc e [l,n], u&eS = Y eAt u^; 

- err = maxj- e n, ||-Une2 - "x^ll- 

We need to store U new , an array of size m x n. 

(3) Only if estimate is . true . : 

Computation of critical At* II: We compute successively: 

- e 2 = maxjgn, \\unil - w 2 ' ^lli 

- Co using (|5 . 10|) with ex and e2; 

- estimate At* oixt of (|5.4|) and set At* = (At* with security factor 
< C < 1 close to one; 

- estimate is set to .false.. 

- If At <£ [/3 At*, 7 At*] with < fi < 7 < 1: estimate is set to 
. true . . 

(4) Only if estimate is . true . and i > 0: 

Computation of s: According to (15.41) with err, Cq and At* = At: 

- e = min{#£, £ mQX } with 9 > 1 as security factor; 

- computation of At* with new e; 

- estimate is set to .false.. 

(5) Computation At new : According to ([2. lip with security factor < 
v < 1 close to one. 

- If At > At*: set err = tol + C with C > 1. Used to potentially 
reject initial At = At . 

- If At™"™ > At* and e 7^ e max : estimate is set to .true.. 
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- At = mm{At new , At*}. 

- If err < tol: t = t + At, i = i + 1, At = min{At,T - t} and 

U — U ne w • 

In this strategy, reaction is always integrated point by point if the reactive term 
is modeled by a system of ODEs without spatial coupling. This integration can 
be performed completely in parallel [9| 7] . On the other hand, for linear diffusion 
problems, another alternative considers a variable by variable resolution, for each 
j€ [l,rn]U^: 

u^=X At u^\ (5.12) 

that can also be performed in parallel [9]. 

Depending on the problem, either the computation of critical At* (steps (1), 
(3) and (4)), or the computation of e (step (4)) can be potentially removed if one 
considers large enough e and fine enough rj. Finally, the whole strategy with all 
steps needs to store at worst two arrays of size I x n and other two of size m x n, 
beyond memory requirements of diffusion and reaction solvers. 

6. Final Numerical Evaluation of the Method 

In this last part, we evaluate the performance of the method in terms of accuracy 
of the simulation, and show that an effective control of the simulation error is per- 
formed in the context of two different problems. First, we will consider a propagating 
wave featuring time-space multi-scale character. Then, the potential of the method 
is fully exploited for a more complex configuration of repetitive gas discharges gen- 
erated by high frequency pulsed applied electric fields followed by long time scale 
relaxation, for which a precise description of discharge and post-discharge phases is 
achieved. 



6.1. BZ equation revisited 

Coming back to BZ model, we perform a time integration of (|5.1ip with several 
accuracy tolerances r\. First of all, we consider the numerical strategy detailed in 
§ 15. 31 without taking into account steps (1), (3) and (4), that is without computation 
of neither critical At* nor e. We set At = 10~ 7 and Eq = 0.05 in all cases, with 
t £ [0, 2]. In this example, a rather small initial splitting time step is chosen to avoid 
initial rejections even though this initial rejection phase usually does not take many 
steps as it will shown in the next example. On the other hand, we have chosen a 
intermediary value for e in order to clearly distinguish the different behaviors of 
the strategy in terms of prediction of the local errors depending on the proposed 
tolerance. 

Figure \5\ shows time evolution of accepted splitting time steps At. In this case, 
BZ equation models a propagating self-similar wave, so splitting time step stabilizes 
once the overall phenomenon is solved within the prescribed tolerance rj. Local error 
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estimates err are also shown, which naturally verify prescribed accuracy, since we 
impose time steps for which err is limited by r/ through (I2.1ip . 




Fig. 5. BZ equation. Time evolution of accepted splitting time steps At (left) and local L 2 error 
estimates err = ||S' At «o — •S'^'^oIIl 2 (right), for several tolerances rj and s = 0.05. 



Table summarizes global L 2 errors between splitting and reference solutions 
at the end of the time domain of study, t — 2. For a fine enough r\ and consequently, 
small enough time steps, a precise error control is achieved by the local error control 
strategy as we could have expected from previous results in Figure [3] for e = 0.05. 
Nevertheless, for 77 = 10 -4 we can see rather high global errors even if this configu- 
ration considers naturally less time integration steps and thus, less accumulation of 
local approximation errors. If we take a look at Figure EJ we note that for s = 0.05 
and local errors of about 10 -4 , the local error estimate, err, is not predicting prop- 
erly real Strang errors, as it was previously discussed, since At > At* . Therefore, 
a strategy that introduces a more precise description of errors for a larger range of 
time steps must be considered, whenever the required accuracy casts the method 
away from its asymptotic behavior. This is an under covered difficulty of any time 
adaptive technique based on a lower order embedded method, and to our knowl- 
edge, an open problem that has not been studied much, and that this work tries to 
overcome. 



Table 2. BZ equation. L 2 errors at final time t = 2 for a, b, c variables and several tolerances rj. 



V 




L 2 error 


a 


L 2 error 


b 


L 2 error 


c 


10" 


-i 


7.97 x 10 


-3 


1.07 x 10 


-2 


4.72 x 10 


-3 


10" 


-6 


1.71 x 10 


-6 


1.83 x 10 


-6 


7.98 x 10 


-7 


10- 


-8 


1.45 x 10 


-8 


1.54 x 10 


-8 


6.78 x 10 


-9 


10- 


10 


1.74 x 10- 


-10 


1.75 x 10" 


-10 


1.08 x 10" 


-10 
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Let us now consider the entire strategy with all steps for several tolerances with 
At = 5x 10 -7 and Eq = 0.05. In the following illustrations we have considered the 
following parameters: e m ax = 0.999; a\ = 1, h\ = ci = a.2 = 1/2, 62 = 2/5 and C2 
1 1 for intermediary time steps evaluations; ( = 0.9 as security factor of critical 
At* estimate; f3 = 0.1 and 7 = 0.95 to define the working region At £ [/3At* , 7 At*] ; 
9 = 10 as security factor of e estimate; C — 10 to potentially reject initial time step 
At ; and v — 0.9 as security factor of At new estimate. All local estimators, err, e\ 
and e2, are computed with normalized L 2 norms. 

Considering the propagating phenomenon, we set N = 10, but we estimate At* 
only twice for i = and i — N. Figure [6] shows time evolution of splitting time 
steps; there are different scenarios depending on the required accuracy. In all cases 
for £0 = 0.05, we estimate initially At* ps 1.4 x 10 -4 . For 77 = 10 -4 , this limitation 
implies smaller time steps than what is required for the prescribed tolerance. Thus, 
At increases until At new > At* and a new e is estimated: e « 0.43. No substantial 
changes are made when i = N, since At £ [/3 Arj* , 7 Af*] for the current r\. 
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Fig. 6. BZ equation. Time evolution of accepted splitting time steps At (left) and local L 2 error 
estimates err = ||S A *«o — S^uoW^ (right), for several tolerances n, considering critical At* and 
computation of e. 



For rj = 10~ 6 , we keep initial At* and £0 since At £ [f3 At* , 7 At*] as we can 
see in Figure [3] Finally, for 77 = 10~ 8 and 77 = 10~ 10 , At < /3At* and thus, e is 
recomputed, giving respectively e ~ 0.016 and 0.0016. In particular, we consider 
larger splitting time steps for which Strang local errors are better predicted. Tabic 
[3] shows that error control is this time guaranteed for all values of tolerance 77, and 
thus, for a larger range of time steps. Compared with previous results in Tabled we 
correct completely the errors in the prediction of local error estimates, which yields 
more accurate resolutions for the largest tolerances; whereas slightly less accurate 
results are obtained for the smallest tolerances since larger splitting time steps are 
considered. 
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Table 3. BZ equation. L 2 errors at final time t = 2 for a, b, c variables and several tolerances r), 
considering critical At* and computation of e. 



V 




L 2 error 
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L 2 error 


b 


L error 


c 


10" 


-4 


6.85 x 10 


-5 


9.04 x 10 


-5 


4.06 x 10 


-5 


10" 


-6 


1.71 x 10 


-6 


1.83 x 10 


-6 


7.98 x 10 


-7 


io- 


-8 


4.53 x 10 


-8 


4.84 x 10 


-8 


2.12 x 10 


-8 


lcr 


10 


4.48 x 10 


-9 


4.77 x 10 


-9 


2.15 x 10 


-9 



6.2. Simulation of multi-pulsed gas discharges 

In this section, we consider a simplified model of plasma discharges at atmospheric 
pressure for which we analyze the performance of the proposed numerical strategy 
in a configuration of nanosecond repetitively pulsed discharges. This kind of phe- 
nomenon is studied for plasma assisted combustion or flow control, for which the 
enhancement of the gas flow chemistry or momentum transfer during typical time 
scales of the flow of 10 -4 — 10~ 3 s, is due to consecutive discharges generated by 
high frequency (in the kHz range) sinusoidal or pulsed applied voltages [23]. As a 
consequence, during the post-discharge phases of the order of tens of microseconds, 
not only time scales are very different from those during discharges of a few tens 
of nanoseconds, but a complete different physics is taking place. Then, to the rapid 
multi-scale configuration during discharges, we have to add other rather slower 
multi-scale phenomena in the post-discharge, such as recombination of charged 
species, heavy-species chemistry, diffusion, gas heating and convection. Therefore, 
it is very challenging to efficiently simulate this kind of highly multi-scale prob- 
lems and to accurately describe the physics of the plasma/flow interaction between 
consecutive discharge/post-discharge phases. 

General model to study gas discharge dynamics is based on the following 
drift-diffusion equations for electrons and ions, coupled with Poisson's equation 

mm- 

d t n c — <9 X • n c v c - <9 X • (D c d x n e ) 
d t n p + <9 X • rtpVp - <9 X • (D p <9 x n p ) 
d t n n - <9 X • ?i n v n - <9 X • (D n <9 x n„) 

e dlV = 

where x G ]R d , n, is the density of species i (e: electrons, p: positive ions, n: negative 
ions), V is the electric potential, Vj = /i^E (E being the electric field) is the drift 
velocity. Di and /x$ , are diffusion coefficient and absolute value of mobility of charged 
species i, q e is the absolute value of electron charge, and eo is permittivity of free 
space, a is the impact ionization coefficient, r] stands for electron attachment on 
neutral molecules, (3 ep and /3 np accounts respectively for electron-positive ion and 



n a|v | - n e ?7|Ve| + n c n p f3 cp + n D j, 
n c a|v c | - n c n p j3 cp + n n n p [3 npi 
n e T]\v c \ - n D n p f3 np - n D -f, 



(6.1) 



-q e (n p 



(6.2) 
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negative-positive ion recombination, and 7 is the detachment coefficient. Electric 
field E and potential V are related by 

E = -d x V. (6.3) 

Nevertheless, in this paper, we will consider a simplified reaction-diffusion ID 
model based on (|6.1[) : 

d t n c - Ddln c = n c a|v c | - n c ?7|v | + n c n p f3 cp , "| 

d t n p - D d%n p = n c a\v c \ - n c n p /3 cp + n n n p (3 np , > (6.4) 

d t n n - D dln D = n ?7|v c | - n D n p f3 np . ) 

As in (|6.ip . all the coefficients of the model are functions of the local reduced 
electric field E/N sas , where E is the electric field magnitude and N gas is the air 
neutral density. Transport parameters and reaction rates for air are taken from 
[2"T] . with attachment coefficients taken from |19j . 

In this numerical illustration, we consider an air gap of 0.5 cm where we have a 
high initial distribution of electrons and ions over the region [0, 0.01] cm. A constant 
electric field of ~ 40 kV/cm is then applied over this region during 10 ns with a pulse 
period of 1 (is. All parameters in (|6.4[) are computed with the imposed field without 
solving neither (16.2[) nor (|6.3p . Finally, we consider a constant diffusion coefficient: 
D = 50cm 2 /s and a spatial discretization of 1001 points. Figure [7] shows the spatial 
distribution of electron density just before and after each pulse. Globally, there 
are at least two completely different physical configurations given either by high 
reactive activity whenever the electric field is applied, or rather by the propagative 
nature of the post-discharge phase. 



0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 

a [cm] x [cm] 

Fig. 7. Repetitive gas discharge model. Spatial distribution of electron density before (left) and 
after (right) each pulse, starting from initial distribution (left) and for a duration of ten pulses. 



Considering the adaptive strategy described in § l5.3l with At = 10 -10 , eo = 0.05 
and the same parameters used for the previous BZ simulation, computation is ini- 
tialized with a time step included in the pulse duration. Figure [5] shows the corre- 
sponding splitting time steps for a tolerance of rj = 10 -3 . Splitting time step features 
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a periodic behavior and succeed to consistently adapt itself to the discharge/post- 
discharge phenomena. This yields high varying time steps going from ~ 10~ 10 to 
~ 10~ 7 . Therefore, after each post-discharge phase, since the new time step is com- 
puted based on the previous one according to (|2.11l) , this new time step will surely 
skip the next pulse. In order to avoid this, each time we get into a new period, 
we initialize time step with the length of the pulse: At = 10ns; this time step is 
obviously rejected as seen in Figure 13 as well as the next ones, until we are able to 
retrieve the right dynamics of the phenomenon for the required accuracy tolerance. 
No other intervention is needed neither for modeling parameters nor for numerical 
solvers in order to automatically adapt time step to describe the several time scales 
of the phenomenon within a prescribed accuracy. 
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Fig. 8. Repetitive gas discharge model. Time evolution of accepted and rejected splitting time 
steps, and imposed electric field for t £ [0, 10] fis (top left), during pulse t £ [5, 5.01] fis (top 
right) and post-discharge t £ [5.01, 6] fis (bottom left). Bottom right: global L 2 errors at the end 
of the pulse (t = 5.01 fis) and the post-discharge phase (t = 6 fis) with and without At* and s 
computation. 



For this application, we compute critical At* and possibly e, for N = 10 and 
TV = 100 in each period in order to perform these computations at least once 
during the discharge and post-discharge regimes. For example, for t £ [5, 6] /zs as 
in Figure [5J e = e max with At* ~ 4.3 x 10~ 9 during the pulse, and e ~ 0.26 
with At* w 1.6 x 10 -7 for the rest of the period. Similar values are found for the 
other periods. Notice that after each pulse, At* is automatically updated because 
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At increases and then At gets equal to At*. In particular, the important difference 
between At* for each region, comes naturally from the completely different modeling 
parameters and hence, physics description of each regime. 

An effective error control is achieved for each part of the phenomenon, as we can 
deduce from the global error between splitting and reference solutions at the end of 
the pulse (t = 5.01 fts) and at the end of the post-discharge phase (t = 6 /xs). If we 
compare these results with the ones obtained without estimating neither At* nor 
e with e = £o = 0.05, we can draw the same conclusions as in the BZ application. 
For less accurate resolutions with high tolerances, the proposed strategy corrects 
the error in the local error estimates made with e = Sq = 0.05; in particular, for 
7] = 10~ 3 there is a ratio of about 10 between both solutions. For higher tolerances, 
-q > 10~ 2 , both methods yield a time step equal to the pulse duration, At = 10ns. 
On the other hand, for the smallest tolerances, slightly more accurate solutions are 
obtained with a fixed e = £o because smaller splitting time steps are used. 

7. Conclusions 

The present work proposes a new resolution strategy for stiff evolutionary PDEs 
based on an efficient splitting scheme previously developed |9|8| that considers high 
order dedicated integration methods for each subproblem in order to properly solve 
the fastest time scales associated with each one of them, and in such a way that 
the main source of error is led by the operator splitting error. Then, to control 
the error of the resolution, it relies on an adaptive splitting time technique that 
allows to discriminate the global time scales related to the coupled phenomenon, 
given a required level of accuracy of computations. Compared with a standard 
procedure for which accuracy is guaranteed by considering time steps of the order 
of the fastest scale, the error control featured by our method implies an effective 
accurate resolution for problems modeling various physical scenarios, independent 
of the fastest physical time scale, and an important improvement of computational 
efficiency whenever highly unsteady phenomena is simulated. In particular, we have 
successfully applied the proposed strategy to a simplified model of plasma discharges 
that nevertheless exhibits a broad time scale spectrum coming from the modeling 
equations and also important and discontinuous variation of parameters in time and 
in space that notably increase the numerical complexity of the problem. 

A numerical analysis of the method has been developed in order to settle a 
solid mathematical background, and a complementary numerical procedure was con- 
ceived in order to overcome classical restrictions of adaptive time stepping schemes 
whenever asymptotic estimates fail to predict the dynamics of the problem. A both 
mathematical and numerical detailed study of the method has thus led to a fully 
complete adaptive time stepping strategy that guarantees an effective control of the 
errors of integration for a large range of time steps; a key issue for problems for 
which splitting time steps can go beyond the fastest physical scales of the problem. 
The contribution of this paper is then mainly given by a dedicated adaptive time 
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splitting method for stiff PDEs, and by a complete study of the behavior of time 
stepping schemes based on lower order embedded methods, for the whole set of 
potential time steps. In this paper we have always considered fine enough spatial 
discretizations in order to perform an evaluation of the theoretical estimates intro- 
duced for the proposed time integration scheme. For higher dimensional problems, 
fine spatial discretization becomes a critical issue in terms of computational costs 
and a technique of local grid refinement might be a good solution to guarantee the 
theoretical behavior of the splitting schemes (see for instance [8]). Nevertheless, a 
mathematical study on the splitting errors with discretized operators will certainly 
be an useful tool to yet improve the performance of these techniques. This and other 
related theoretical aspects are particular topics of our current research. 
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